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Abstract 

We give the first exact algorithmic study of facility location problems that deal with finding 
a median for a continuum of demand points. In particular, we consider versions of the "contin- 
uous fc-median (Fermat- Weber) problem" where the goal is to select one or more center points 
that minimize the average distance to a set of points in a demand region. In such problems, 
the average is computed as an integral over the relevant region, versus the usual discrete sum of 
distances. The resulting facility location problems are inherently geometric, requiring analysis 
techniques of computational geometry. We provide polynomial-time algorithms for various ver- 
sions of the L\ 1-median (Fermat- Weber) problem. We also consider the multiple-center version 
of the Li fc-median problem, which we prove is NP-hard for large k. 
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1 Introduction 



" There are three important factors that determine the value of real estate - location, 
location, and location." 

The Fermat- Weber Problem. There has been considerable study of facility location problems 
in the field of combinatorial optimization. In general, the input to these problems includes a 
weighted set D of demand locations, with weight distribution 8 and total weight fi, a set F of 
feasible facility locations, and a distance function d that measures cost between a pair of locations. 
In one important class of questions, the problem is to determine one or more feasible median 
locations c G C C F in order to minimize the average cost from the demand locations, p £ D, to 
the corresponding central points c p E C that are nearest to p: 



where d(p, C) = mxn c ^c d(p, c). If there is one median point to be placed, the problem is known as 
the classical Fermat-Weber problem; its history reaches back to Fermat, who first posed it for three 
points, a case that was first solved by Torricelli. (Note that this special case has another natural 
generalization: in the well-known Steiner tree problem, the objective is to find a connected network 
of minimum total length connecting a given set of points. See [23] for a recent study of the relation 
between these problems and further discussion.) In the context of facility location, the median 
problem was discussed in Weber's 1909 book on the pure theory of location for industries [HH] (see 
|61| for a modern survey); because of this connection, we will speak of the Fermat-Weber problem 
(FWP) throughout this paper. More generally, for a given number k > 1 of facilities, the problem 
is known as the k-median problem. A problem of similar type with a different objective function is 
the so-called k-center problem, where the goal is to find a set of k center locations such that the 
maximum distance of the demand set from the nearest center location is minimized. 

Geometric Facility Location. There is a vast literature on location theory; for a survey, see 
the book of Drezner |18j . with its over 1200 citations that not only include papers dealing with 
mathematical aspects of optimization and algorithms, but also various applications and heuristics. 
A good overview of research with a mathematical programming perspective is given in the book of 
Mirchandani and Francis |45j . 

With many practical motivations, geometric instances of facility location problems have at- 
tracted a major portion of the research to date. In these instances, the sets D of demand locations 
and F of feasible placements are modeled as points in some geometric space, typically !R 2 , with 
distances measured according to the Euclidean (L2) or Manhattan (L\) metric. In these geometric 
scenarios, it is natural to consider not only finite (discrete) sets F of feasible locations, but also 
(continuous) sets having positive area. For the classical Fermat-Weber problem, the set F is the 
entire plane 5i 2 , while D is some finite set of demand points. 

There has been considerable activity in the computational geometry community on facility 
location problems that involve computing geometric "centers" and medians of various types. The 
problem of determining a 1-center, i.e., a point c to minimize the maximum distance from c to 
a discrete set D of points, is the familiar minimum enclosing disk problem, which has linear-time 
algorithms based, e.g., on the methods of Megiddo. The geodesic 1-center of simple polygons has an 
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0(n log n) algorithm 53 j; in this version of the center problem, distances are measured according 
to shortest paths (geodesies) within a simple polygon. Recent results of Sharir et al. [HI I22| 155] 
have yielded nearly-linear-time algorithms for the planar two-center problem. The more general 
p-center problem has been studied recently by [56] . 

Continuous Location Problems. Location theory distinguishes between discrete and contin- 
uous location theory (see |27j). However, for median problems, this distinction has mostly been 
applied to the set of feasible placements, distinguishing between discrete and continuous sets F. 
It is remarkable that, so far, continuous location theory of median problems has almost entirely 
treated discrete demand sets D [23 E2j. We should note that there are several studies in the 
literature that deal with fc-center problems with continuous demand, e.g, see |43l I58j . where de- 
mand arises from the continuous point sets along the edges in a graph. See for results on 
the placement of k capacitated facilities serving a continuous demand on a one-dimensional in- 
terval. Also, fc-center problems have been studied extensively in a geometric setting, see e.g. 
[IllIZll2HllSnilS21ISSllSlllSSllSnill21llllE3Eni- However, designing discrete algorithms for fc-center 
problems can generally be expected to be more immediate than for fc-median problems, because 
the set of demand points that determine a critical center location will usually form just a finite set 
of d + 1 points in o!-dimensional space. 

Continuous demand for fc-median problems is also missing from the classification in We 
contend that the practical and geometric motivations of the problem make it very natural to consider 
exact algorithms for dealing with a continuous demand distribution for /c-median problems: if a 
demand occurs at some position p £ D, according to some given probability density 5{p), then 
we may be interested in minimizing the expected distance f peD d(c,p)5(p)dp for a feasible center 
location c S F. 

To the best of our knowledge, there are only few references that discuss /c-median problems 
with continuous demand: See the papers |51| 165] for a discussion of continuous demand that arises 
probabilistically by considering a discrete demand in an unbounded environment with a large 
number of demand points, leading to a heuristic for optimal placement of many center points. 
Drezner |19j describes in Chapter 2 of his book that normally a continuous demand is replaced 
by a discrete one, for which the error is "quite pronounced for some problems". (See his chapter 
for some discussion of the resulting error.) Wesolowsky and Love |62| (and also in their book |41| 
with Morris) and Drezner and Wesolowsky [20] consider the problem of continuous demand for 
rectilinear distances. Practical motivations include the modeling of postal districts and facility 
design. They compute the optimal solution for one specific example, but fail to give a general 
algorithm. More recently, Carrizosa, Muhoz-Marquez, and Puerto [Tj |H] use convexity properties 
for problems of this type to deal with the error resulting from nonlinear numerical methods for 
approximating solutions. It should be noted that the objective function is no longer convex when 
distances are computed in the presence of obstacles. 

In this paper, we study the fc-median problem, and its specialization to the Fermat- Weber 
problem (k = 1), in the case of continuous demand sets. Another way to state our continuous 
Fermat-Weber (1-median) problem is as follows: In a geometric domain (e.g., cluttered with ob- 
stacles), determine the ideal "meeting point" c* that minimizes the average time that it takes an 
individual, initially located at a random point in D, to reach c*. Another application comes from 
the problem of locating a fire station in order to minimize the average distance to points in a neigh- 
borhood, where we consider the potential emergencies (demands) to occur at points in a continuum 
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(the region defining the neighborhood D). As we noted above, this objective function is different 
from the situation in which we want to minimize the maximum distance instead, a problem that 
has been studied extensively in the context of discrete algorithms. 

Choice of Metric. Many papers on geometric location theory have dealt with continuous sets 
F of feasible placements, including [21 El d31 [T^l [THl Ell EZl EH1 EH1 EH1 - In the majority of 
these papers, distances are measured according to the L\ metric. In fact, it was shown by Bajaj [S] 
that if L>2 distances are used, then even in the case of only five demand locations (\D\ = 5), the 
problem cannot be solved using radicals; in particular, it cannot be solved by exact algorithmic 
methods that use only ruler and compass. (Chandrasekaran and Tamir ^Jll give a polynomial-time 
approximation scheme that uses the ellipsoid method.) In this paper, we too concentrate on the 
problem using the L\ metric. While we can exactly solve some very simple special cases in the L2 
metric, in general the integrations that are required to solve the problem are likely to be just as 
intractable as the classic Fermat- Weber problem. 

Summary of Results. In this paper, we give the first exact algorithmic results for location 
problems that are continuous on both counts, in the set D as well as the set F. In our model D 
and F are each given by polygonal domains. Our goal is to compute a set of k (k > 1) optimal 
centers in the feasible set F that minimize the average distance from a demand point of D to the 
nearest center point. Our results include: 

(1) A linear-time (O(n)) algorithm for computing an optimal solution to the 1-median (Fermat- 

Weber) problem when D = F = P, a simple polygon having n vertices, and distance is taken 
to be L\ geodesic distance inside P. 

(2) An 0(n 2 ) algorithm for computing an optimal 1-median for the case that D = F = P, a 

polygon with holes, and distance is taken to be (straight-line) L\ distance. 

(3) An 0(1 + nlogn) algorithm (where I = 0(n 4 ) is the complexity of a certain arrangement) 

for computing an optimal 1-median for the case that D = F = P, a polygon with holes, and 
distance is taken to be L\ geodesic distance inside P. 

(4) A proof of NP-hardness for the £>median problem when the number of centers, k, is part of the 

input, and D = F = P is a polygon with holes. This adds specific meaning to the statement 
by Wesolowsky and Love that computing the optimal position of several locations "is 
obviously very tedious when (the number of locations) is very large" . 

(5) Generalizations of our results to the following cases: non-uniform probability densities over 

the demand set D; fixed-orientation metrics (generalization of L%), which can be used to 
approximate the Euclidean metric; higher dimensions; and, F 7^ D for straight-line distances. 

This paper represents research done as part of the PhD thesis of Weinbrecht; additional exam- 
ples, discussion, and details may be found in [HO] . 

2 Preliminaries 

Basic Definitions. We will let Z = (x,y) denote a candidate center point in F C 3ft 2 . (We 
concentrate on two-dimensional problems until Section |HJ where we discuss extensions to higher 
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Figure 1: A simple example in which the region D = F = P is a polygon (square) with a single 
(parallelogram) hole, shown shaded. The average straight-line L\ distance is minimized by two 
optimal center points, each marked "c*", on the boundary of P. This example is analyzed in more 
detail later; refer to Figure HUP 



dimensions.) We defer discussion of multiple center points {k > 1) to Section [7J for now, k = 1 and 
we consider the Fermat- Weber (1-median) problem. 

We let P denote a polygonal domain: This is a connected planar set of points that is bounded by 
a finite set of disjoint simple closed polygonal curves. We assume that P is nondegenerate, i.e., it is 
a closed set that equals the closure of its interior points; in particular, the interior is connected. We 
say that the vertices of P are the vertices of its boundary; as part of the nondegeneracy assumption, 
we assume that each vertex is incident to precisely one boundary polygon and to two edges. In the 
case of one connected boundary, we say that a polygon P is simple; otherwise we say that P has 
one or several holes, i.e, bounded components of the complement 5ft 2 \ P. A critical vertex v of P is 
one that has locally extremal x- or y-coordinate relative to the boundary component containing v, 
and an interior angle of at least tt. A chord of P is a straight line segment within P that connects 
two points on the boundary of P. If P is a simple polygon, any chord subdivides P into two or 
more pieces. 

For purposes of our discussions, we focus on the case in which D = F = P: we restrict Z to 
P, which also equals the demand set. Furthermore, we focus our discussion to the case in which 
the demand is uniformly distributed over the set D = P, so our goal is to minimize the average 
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distance, f(Z), given by the integral 

f{Z) = f(x,y) = -[ [ d(Z,(u,v))dudv, 
J (u,v)eP 

where fi is the total area of P, d{-, •) denotes either (straight-line) L\ distance or geodesic (shortest- 
path) L\ distance within P. (We abuse notation slightly by writing f(Z) = f((x,y)) = f(x,y).) 

Shortest Path Maps. In order to analyze the fe-median problem with respect to geodesic dis- 
tances, we will utilize several definitions and results from the theory of geometric shortest paths 
among obstacles; see Mitchell [131301 for surveys on the subject of geometric shortest paths. 




Figure 2: An example of the L\ geodesic SPM(Z) for a polygon P with two holes (shaded). Vertices 
are labeled with their L\ geodesic distances from Z . The cell rooted at Z is partitioned into four 
quadrants by a vertical and a horizontal chord through Z. Within each other cell, rooted at a vertex 
of P, an arrow is drawn pointing to the root of the cell. Vertices of P whose cells are nonempty 
are drawn as solid black dots; the white vertices have empty cells in SPM(Z). 

For a given geometric environment P, the shortest-path map, SPM(Z) with respect to source 
point Z represents the set of shortest paths within P from Z to all other points t £ P. The SPM(Z) 
is a decomposition of P into cells a, each having a unique root vertex, r a , such that a shortest path 
within P from Z to t £ P is given by following a shortest path within P from Z to r a , then going 
from r a to t directly along a "straight" segment. A straight segment between p £ P and q £ P 
is a path whose length is d(p,q), where d(-,-) denotes our underlying distance function (L\, L2, 
etc.). If the underlying distance is L2, a straight segment is necessarily a straight line segment; 
if the underlying distance is L\, as in most of this paper, a straight segment from p to q consists 
of any path from p to q that is both x-monotone and y-monotone (e.g., an "L-shaped" path or a 
"staircase" is straight in the L\ metric). If t lies in a cell a rooted at r a , the geodesic distance from 
Z to t is given by dc(Z, t) = dc(Z, r CT ) +d(r a ,t), where da(-, •) denotes the shortest-path (geodesic) 
distance function induced by P. There is also a cell with root Z, consisting of points t £ P for 
which a shortest path from Z to t is attained by a straight segment Zt. Points t on the common 
boundary of two or more cells are optimally reached by shortest paths to t whose last segment 
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joins t to the root vertex of any one of the cells having t on its boundary. Refer to Figure |2] for 
an example of SPM(Z) in the L\ geodesic metric. Each vertex v of P is the root of a (possibly 
empty) cell. A vertex v = r a that is the root of a cell a lies on the boundary of a neighboring 
cell a' , so there is a shortest path from Z to r a whose last segment is r a /r a ; we say that r a < is 
the predecessor of r a . In addition to its predecessor, we store with each vertex v of P the length, 
da(Z,v), of a shortest path within P from Z to v. In summary, SPM(Z) is a decomposition of P 
into cells according to the "combinatorial structure" (sequence of obstacle vertices along the path) 
of shortest paths from Z to points in the cells. 

For any two distinct vertices u and v of P, the bisector with respect to u and v is the locus of 
points p G P that are (geodesically) equidistant from Z via the two distinct roots, u and v; thus, 
a bisector is the (possibly empty) locus of points p £ P satisfying dc{Z,p) = dc(Z,u) + d(u,p) = 
da(Z,v) + d(v,p). If the underlying metric is Euclidean, bisectors are curves, easily seen to be 
straight line segments or hyperbolic arcs. In the L\ geodesic metric, bisectors may be horizontal 
line segments, vertical line segments, polygonal chains of segments that are horizontal, vertical, or 
diagonal (with slope ±1), or they may, in a degenerate situation, consist of regions. See Figure El 
and Figure 0] In order that cells of the SPM(Z) do not overlap in regions of nonzero area, and so 
that the SPM(Z) is a planar decomposition of P, it is convenient to resolve degeneracies so that 
all bisectors are one-dimensional (polygonal) curves, not regions. We do this as follows. First, if 
p G P satisfies do(Z,p) = dc{Z^u) + d(u,p) = da(Z,v) + d(v,p), then we consider p to be in 
the cell rooted at u (and not in the cell rooted at v) if d(u,p) < d(v,p); a consequence is that 
the points p£P that are in the cell rooted at u are visible to u (i.e., the segment up lies in the 
cell). Figure shows the result of applying this rule to two situations in which the L\ geodesic 
bisector would otherwise be a region. Second, we innnitesimally perturb the vertices of P so that 
no two vertices lie on a common line of slope ±1; this implies that set of points p G P for which 
da(Z,p) = dc(Z,u) + d(u,p) = dc(Z,v) + d(v,p) and d(u,p) = d(v,p) is a polygonal chain, not 
a region of nonzero area. (All of our results apply also to the unperturbed problem instances, by 
standard arguments.) 

It is easily seen that a small horizontal (or vertical) shift of Z by an amount h results in a shift 
of the bisector between two vertices u and v by an amount h/2. 




(a) (b) (c) 



Figure 3: Examples of the L\ bisectors between two points, Pi and P2, in the L\ metric. In each 
of the three cases, the locus is shown of points p for which d(P\,p) = d(P2,p). Note that in case 
(c), when Pi and P2 lie on a line of slope ±1, the bisector includes two regions (shown shaded). 

Some bisector curves, such as those horizontal and vertical segments shown solid in Figure [2J 
may be crossed by shortest paths from Z to points p G P. However, bisector curves that consist of 
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fa) fbl 

Figure 4: Examples of L\ geodesic bisectors. In both cases, the points in the lightly shaded quadrant 
are equidistant from Z via u and via v; our tie-breaking rule, however, assigns the shaded bisector 
points to the cell rooted at u, resulting in the bisector curve shown dashed. 



points that are the endpoints of maximal shortest (geodesic) paths are not crossed by any shortest 
path; these are shown dashed in Figure |2] and are called watersheds. (A shortest path from Z to 
p € P is maximal if it is not the proper subset of a shortest path from Z to some other point p' ^ p 
of P.) One can think of the watershed bisectors as "ridges" that partition P into regions according 
to the topological type (homotopy class) of shortest paths; a point on a watershed can be reached 
from at least two homotopically distinct shortest paths from Z. 

Shortest path maps can be computed in optimal time 0(n log n) for a polygon with n vertices, 
both in the Euclidean metric and in the L\ metric |29| I46| 13%]. More generally, shortest-path maps 
can be applied to weighted region metrics, where the time for traveling depends on a local density 
function. For more information, see the surveys of Mitchell |49| lf>0] . 

Finally, we define one other piece of notation. For any position Z of a center, we define two 
subsets of the domain P: the set W{Z) (resp., E{Z)) of all points p G P for which every shortest 
path from Z to p enters the open halfplane to the left (resp., right) of Z prior to entering the open 
halfplane to the right (resp., left) of Z. In other words, W(Z) (resp., E(Z)) corresponds to the set 
of points p £ P for which an optimal path to p initially heads to the west (resp., east). Similarly, 
P is subdivided into the sets N(Z) and S(Z) of points for which shortest paths initially head 
north versus south. These definitions apply both to geodesic distances, and to straight- 

line distances; refer to Figure El for illustrations. While there may be points p G P that are not 
in W{Z) U E(Z), such points either lie on the vertical line through Z or, in the case of geodesic 
distances, lie on a bisector (since such points are reached by at least two distinct optimal paths - one 
heading initially east, one heading initially west). By the convention and perturbation assumption 
that allows us to assume that bisectors are one-dimensional curves (rather than regions), we see 
that W(Z) U E(Z) covers all of P except for a set of zero area; a similar statement holds for 
N(Z) U S{Z). In the following, we will use w(Z), e(Z), n(Z), s(Z) for the area of W(Z), E(Z), 
N(Z), S(Z), respectively. 

3 Local Optimality Conditions 

For any given center location Z, the objective function value f{Z) that gives the average distance 
from Z to points of P can be evaluated by decomposing P into a set of "elementary" pieces, 
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computing the average distance for each piece, and then obtaining the total average distance as a 
weighted sum of the average distances for the pieces. In the case of straight-line L\ distance, we 
simply use a trapezoidization (or triangulation) of P to determine the pieces; this can be done in 
linear time if P is simple, or in 0(n log n) time if P has holes. In the case of geodesic distance, the 
shortest-path map, SPM(Z), gives a decomposition of P into cells (each of which can be refined into 
triangles or trapezoids to yield a decomposition into 0(l)-size pieces), each having a corresponding 
root vertex on its boundary. By computing the average distance from points of a cell to the cell's 
root r and adding this average to the distance dc(Z, r), and then summing over all cells, we obtain 
the average geodesic distance, f(Z). 

The average distance associated with a single elementary piece is given by the following result, 
which can be verified easily by straightforward integration. As any region can be subdivided into 
a limited number of triangles of this type, it can be used as a stepping stone for computing the 
objective value for more complicated regions. 



C 



y 




Figure 5: LemmaUnotation used in computing the average L\ distance of AABC from the point A. 



Lemma 1 For a triangle r with vertices A, B, and C , such that edge AB is horizontal, let a be the 
length of AB, let c be the length of the altitude from C, and let b be the distance from A to the foot 
of the altitude from C. Then the average L\ distance of points in r from vertex A is 3 (a + b + c) . 

The objective function, f(Z), is a continuous function of the location of the center Z. Because 
the set F = P of feasible placements is a compact domain, it follows that there is an optimum. If Z 
is a point minimizing f(Z), then, Z must be locally optimal, meaning that there cannot be a feasible 
direction H = (xh,Uh), i-e., Z + eH G P for sufficiently small e, such that f(Z) > f(Z + eH). 
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If / is differentiable at some point Z £ P, then (Vf(Z),H) > 0. In particular, for interior points 
that are locally optimal, the gradient must be zero; for points in the interior of boundary edges, 
the gradient must be orthogonal to the boundary. In the following lemma we compute the gradient 
of/: 

Lemma 2 Consider the objective function f for average straight-line L\ distance in a region P 
of area /i = /x(-P). Let Z = (x,y) be a point in P. Then the first partial derivatives of f are 
well-defined and given by: 

f x (Z) = ±(w(Z)-e(Z)), 

f y (Z) = Us{Z)-n{Z)). {L) 




Proof. We compute f x ; f y is computed similarly. Refer to Figure (left). Consider the point 
Z' = Z + (h, 0) for some sufficiently small h. Let C(Z, Z') := {p = (x p , y p ) G P\x < x p < x + h} be 
the narrow vertical strip between Z and Z'. We compute: 



f(x y) - lim f{x + h >ti ~ fe^ - lim /( ^ = 



( i 



lim 



i / dip, z') dp- J d(p,z)d P 



v P eP 



h 
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(\[ J d(p,Z')dp- J d(p,z)d P ^ Xx 
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oew{z) 



peW(Z) 
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(\[ J d(p,Z')dp- J d(p,Z)dp^ S 
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+ lim 



oec{z,z') 



p£C(Z,Z') 
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V 

(\ J d(p,Z')dp- J d(p,Z)dp^ 



pee(z') 



p&E(Z>) 



\ 



h 



For small enough h, the area of C(Z,Z') is a linear function of h, since the boundary of P 
is made up of straight line segments. The term involving the difference of the two integrals over 
the points p E C(Z,Z') is dependent only on the x-contribution to the L\ distance (\x p — xz\ or 
\x p — %z'\) between p and Z or Z', since the y-contributions cancel (yz = Uz 1 )- Since x p varies 
within a range of h, and \x p — xz\ < h (and \x p — xz>\ < h) as well, we get that 



\pec(z,z 



') P ec(z,z') / 



Since 



and 



we get 



d(p,Z')dp= J (h + d(p,Z))dp, 
pew(z) pew(z) 

J d(p,Z)dp= J {h + d(p,Z'))dp, 

p£E{Z) p£E(Z) 



fx{x z ,y z ) = hm 

h— >0 



/ / hdp- J hdp + 0(h 2 ) 
| pew(z) 



peE(Z) 



-w(Z)-e(Z) 



as claimed. 



□ 



The above lemma characterizes the gradients in the case of straight-line distances. For geodesic 
distances, we proceed in a similar manner; refer to Figure El (right). Note that the shape of C(Z, Z') 
is different in the case of geodesic distances. For the example in the figure, C(Z, Z') is the union of 
one connected strip (shown white), bounded by vertical line segments through Z and Z' and the 
boundary of P, and possibly several narrow regions that are swept by the watersheds as Z moves 
to Z' (shown darkly shaded in the figure); note that these latter regions only occur in the case of 
P being a polygon with holes. 
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In the lemma below, we state the result under a technical assumption, which avoids difficulties 
with the continuity of w(Z) and e(Z), as well as with degenerate bisectors. 

Lemma 3 Consider the objective function f for average geodesic L\ distance in a region P of 
area /x. Let Z = (xz,yz) be a point in P; assume that neither xz nor yz coincide with the x- or 
y-coordinate of a critical vertex of P, and that Z does not lie on a watershed bisector in a shortest 
path map, SPM(v), with respect to a vertex v of P. Then the first partial derivatives of f are 
well-defined and given by: 

f x (Z) = ±(w(Z)-e(Z)), 

f y (Z) = i (*(Z) - *(Z)) . [Z) 

Proof. The argument is directly analogous to the proof of Lemma |2J The assumption that Z does 
not lie on a watershed bisector of SPM(u) for any vertex v implies that no watershed bisector of 
SPM(Z) contains a vertex v of P. This implies that the area of C(Z, Z') is a linear function of the 
perturbation parameter h. Furthermore, the difference d(p, Z) — d(p, Z') for points p £ C(Z, Z') is 
bounded by h. The claim follows, as in the previous lemma. □ 

In some situations, we make use of properties of higher-order derivatives of /. In particular, we 
use the following lemma: 

Lemma 4 The objective function f is piecewise the sum of two cubic functions, fi(x) and f2(y), 
both for straight-line and for geodesic distances. 

Proof. We start by considering the objective function / for average straight-line L\ distances in P. 
Let Z = (xz, yz) be a point in P, with neither xz nor yz coinciding with the x- or y-coordinate of 
a critical vertex of P. Consider a small change in the x-coordinate of Z: let Z' = Z + (h, 0). Then 
W(Z') = W{Z) U C(Z, Z') and E(Z') = E(Z) U C(Z, Z'), where C(Z, Z') = B U B x U . . . U B k 
is a union of a set of trapezoids P>i, each of width h; refer to Figure El Since the area of each 
trapezoid B>i is a quadratic function of the width h, for small enough h, we get that w(Z') — w(Z) 
and e(Z') — e(Z) are also quadratic in h. Since f x (Z) = j^(w(Z) — e(Z)), by LemmaEl this implies 
that f xxx (Z) is a constant. (Specifically, f xxx (Z) = (2/fi) Yli( m i,t ~ m i,b), where m i)t (resp., m i)6 ) 
is the slope of the edge of P bounding the top (resp., bottom) of trapezoid B>i.) To see that / does 
not contain any terms that have xy as a factor, note that for fixed x, W{Z) and E{Z), and thus 
f x (Z) does not depend on y. 

For geodesic distances, the claim follows in a similar manner. Now, the regions P>i forming the 
connected components of C(Z,Z') are not necessarily vertical- walled trapezoids, but have paral- 
lel walls formed by translates of watershed bisectors. See Figure H3 Instead of being trapezoids 
bounded by vertical chords through Z and Z 1 , however, the areas E>i are bounded by two bisectors, 
corresponding to translates of watershed bisectors for Z and Z' . Specifically, it follows from prop- 
erties of L\ bisectors described in Section |^1 that these bisectors move in a parallel fashion, provided 
that no degeneracy of a bisector occurs during the move from Z to Z 1 ', i.e., no polygon vertex is 
hit by a bisector. The area of region B{ is thus quadratic in h, and the result follows. Again, for 
fixed x, f x (Z) does not depend on y. □ 
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4 Straight-Line L\ Distance 



For straight-line distances, a finite average is guaranteed even for disconnected regions P, as long as 
they are compact. In the following, we will consider local optimality for finding a global optimum 
of /. Lemma 13 motivates considering the L\ origin of P, which is a point Z with w{Z) = e{Z) 
and n(Z) = s(Z), i.e., the (unique) point that is both a median of the x- and the y-distribution. 
However, the example in Figure shows that even for the special case of a simple polygon P, the 
L\ origin of P may not be a feasible point. 







p 










1 













Figure 7: For straight-line distances, the L\ origin of a simple polygon P may be infeasible. 

This makes it slightly involved to compute all local optima. In the following, we describe how 
to evaluate them in 0(n 2 ) time. 

Theorem 5 For straight-line L\ distances, a point Z* = (x*,y*) in a polygonal region P that 
minimizes the average distance f to all points in P can be found in time 0(n 2 ). 

Proof. We apply the local optimality conditions. We start by computing in time O(nlogn) the 
L\ origin Z m of P; if Z m is feasible (i.e., Z m 6 P), we are done. If no interior point of P is a local 
optimum, then we have to consider the boundary of P. This yields a set E of 0(n) line segments 
and a set V of O(n) vertices that we examine for local optimality. 

We overlay the set of vertical and horizontal lines through all vertices of P with E, subdividing 
each segment in E into O(n) pieces, bounded by a total of 0(n 2 ) "overlay" vertices V a . Let E be 
the resulting set of 0(n 2 ) subsegments. See Figure |S1 

Now we can examine the interior points pt = (t,y(t)) of each edge ej E E for local optimality. 
Let Sj be a vector parallel to ej. By construction of E D , the vertical and horizontal lines through 
Pt cannot encounter a vertex of P as pt slides along one subsegment of E Q . Thus, it follows from 
Lemma|I]that f x (pt) is a quadratic function in t, and so is f y (pt)- Therefore, considering for Z £ ej 
the local optimality condition 

<v/(ptW = o 
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for all 0(n 2 ) subsegments ej S E D yields a set of 0(n 2 ) quadratic equations in t. These can be 
solved in amortized time 0(n 2 ), because we can obtain the coefficients of each quadratic equation 
in amortized constant time by advancing from cell to cell in the overlay arrangement. This gives, 
for each subsegment ej, at most two local optima, qj t \ and q^i. Let Vi be the union of E Q and all 
qj t i and q^2- By construction, Vi contains 0(n 2 ) elements, and all local optima of / occur at points 
of V(. Thus, our goal is to evaluate the objective function at each of these points and to select the 
best one. This is simply done in amortized time Oil) per candidate, by walking over the overlay 
arrangement and incrementally updating the value of the objective function. □ 

In many cases, the following property of straight-line medians can be applied for a reduction 
of the set of boundary segments that we need to consider. (See Figure EH for an illustration.) If 
Z m = (x m ,y m ) is the L\ origin of P and p\ = {x\,y\) and P2 = (%2,y2) are points in P, we say 
that pi dominates P2 , if pi lies in the rectangle spanned by Z m and p2 ■ 

y a 




Figure 9: In the region P (shown shaded), P2 is dominated by p± and cannot be a local optimum. 
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Lemma 6 Let pi = (xi,yi) and p2 = (2:2,2/2) be points in P. If pi dominates P2, then f(p\) < 

f{P2). 

Proof. Suppose that %i > x\ > x m and y% > y\ > y m . Then w(jp2) > w{pi) > ^ and 
s(p2) > s (pi) > 2 ' so mov i n § a center from P2 to p\ cannot increase the objective value. □ 

Using a plane-sweep algorithm, it is possible to identify the non-dominated portions of the 
boundary in time O(nlogn). If this set has complexity o(n), then we get a reduction of the overall 
complexity. 

A simple example of the problem solved in this section is shown in Figure EH (based on the 
example given in Figure ^) . This example shows the necessity of solving quadratic equations in 
computing the optimal solutions (points pa and ps). The non-dominated points are pi, P2, and ps, 
as well as the points interior to the segment e2,3 = P2P3- 
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Figure 10: The example from Figure ^ is analyzed for the case of straight-line L\ distances. There 
are two optimal center points: P4 = (2y7 — |, 11 — 4\/7) and its mirror image, ^5. 



5 Geodesic L\ Distances in Simple Polygons 

In this section, we show how to compute in optimal (O(n)) time a point that minimizes the average 
geodesic L\ distance for a simple polygon P without holes. 

From Lemma 01 we know that the partial derivative f x {Z) = j^(w(Z) — e(Z)) vanishes if 
w{Z) = e(Z). Therefore we consider the functions w(Z) and e(Z) that are well-defined even for 
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points for which the gradient is not. As x increases, w(Z) increases monotonically, while e(Z) 
decreases monotonically. Note that w(Z) may be discontinuous at critical vertices of P: As shown 
in Figure an entire region may switch from "east" to "west" as a vertical chord passes through 
a critical vertex. 

However, even discontinuous behavior at critical coordinates does not impair monotonicity of 
w{Z) and e(Z), so there is still a well-defined vertical median chord c x at some x-coordinate x m 
such that w(Zi) < e(Z\) for all Z\ = (x%,yi) with x\ < x m (implying f x {Zi) < just left of 
x m ), and w{Z\) > e{Z\) for all Z2 = (2:2,2/2) with X2 > x m (implying f x {Z?) > just right of 
x m ). Similarly, there is a unique horizontal median chord c y at y-coordinate y m . Again we call 
Z m = (x m ,y m ) the L\ origin of P. 

In the following, we use the structure of simple polygons to show that the locally optimal point 
Z m has to belong to the feasible region P (possibly on the boundary of P), implying that it is a 
unique global optimum. 

Theorem 7 The point Z m is feasible (lies in P) and thus a unique global optimum, minimizing 
the average L\ geodesic distance to points in P. 

Proof. The chord c x subdivides P into two pieces: let E denote the part to the right ("east") 
of c x , and W the part to the left ("west") of c x . Note that E or W may consist of two or more 
connected components only if c x passes through a critical vertex. Similarly, the chord c y subdivides 
P into the region N ("north") that lies above c y , and the region S ("south") that lies below c y . 

We claim that c x and c y intersect at a point (Z m ) inside P. The proof is by contradiction; 
assume that Z m lies outside P. We will distinguish the following cases. 

Case 0: Neither c x nor c y are critical. If Z m P, then simplicity of P implies that the 
two chords subdivide P into three pieces; this means that precisely one of the pieces W and E has 
nonempty intersection with one of the pieces iV and S. Without loss of generality, assume that the 
two chords subdivide P into the three pieces, E, N n W, and S, as shown in Figure ITT1 Since P 
is a nondegenerate polygon, the corresponding areas, (J,(E), /j,(N C\ W), and fJ,(S), are all positive, 
with (jl{E) + n(S) + fi(N n W) = fi(P) = \i. However, the local optimality of x m implies that 
n(E) = ^ and the local optimality of y m implies that //(£) = ^, implying the contradiction that 

/i(JVnw) = o. 

Case 1: Exactly one of c x and c y is critical. Without loss of generality, assume that c y 
passes through a critical vertex of P that is a local maximum of the boundary of P, as in Figure IT21 
As in Case 0, the (noncritical) vertical chord c x partitions P into two pieces, E and W, each of area 
77. Also, Cy partitions P into one "upper" piece N and two "lower" pieces, Si and S2, each with 
positive area. (In degenerate situations, c y may pass through multiple critical vertices, resulting in 
multiple lower pieces; our arguments includes this case by considering additional pieces as part of 
5*2.) The assumption that c x and c y do not cross inside P implies that either W or E is a strict 
subset of Si, S2, or N. This is a contradiction, since fJ>(W) = fJ-(E) = ^, and (i(N) = ^, /x(£i) < f , 
and //(S2) < 2- 

Case 2: Both c x and c y are critical. Without loss of generality, assume that c x passes 
through a critical vertex of P that is a local minimum of the boundary, while c y passes through a 
critical vertex that is a local maximum of the boundary, as shown in Figure As in Case 1, the 
horizontal chord c y subdivides P into a "northern" piece N, and two "southern" pieces, Si and S2, 
each with positive area. Similarly, the vertical chord c x subdivides P into a "western" piece W, 
and two "eastern" pieces, Ei and Ei, each with positive area. We assume further, without loss of 
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Figure 11: Case 0: Neither of the median chords c x and c y is critical. 




Figure 12: Case 1: Exactly one of the median chords is critical. 



generality, that the critical vertex through which c x passes lies within the piece S\, as shown in the 
figure. We have a contradiction in the fact that fJ>{N) = fJ>(W) = ^, while the set P \ (N U W) has 
positive area greater than //(S^)- □ 



Theorem 8 The point Z m can be computed in linear time. 



Proof. We describe how to compute the x-coordinate x m of Z m ; the y-coordinate is found in a 
similar manner. 

In linear time (using Chazelle's algorithm jllj). we build the vertical trapezoidization of P, 
which is defined by drawing vertical chords through every vertex of P. Each piece, r«, of the 
resulting subdivision is either a vertical-walled trapezoid or a triangle having one side vertical 
(such a triangle can be considered to be a degenerate vertical- walled trapezoid). Consider the 
adjacency graph Q of these pieces (i.e., the planar dual of the trapezoidization); because P is a 
simple polygon, Q is a tree. 

Let T m denote the trapezoid containing the vertical chord c x through Z m . (We assume, without 
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Figure 13: Case 2: Both of the median chords are critical. 



loss of generality, that c x is not one of the vertical walls of r m ; the degenerate case is readily handled 
by similar arguments.) Let C max be a connected component of P \ r m that has maximum area; let 
r max be the unique trapezoid within C max that is (vertical wall) adjacent to r m . The area of C ma x 
cannot exceed ///2, by the local optimality criterion. (Moving Z m from r m by an infinitesimal e into 
7"max would reduce the distance to Z m by e for a set of points of a total area more than n/2, while 
increasing it by at most e for a set of points of total area less than fi/2.) Thus, r m corresponds to 
what is called a median node in the weighted tree Q, whose nodes are weighted by the areas of the 
corresponding trapezoids. See Figure ITU 

A median in a weighted tree can be computed in linear time (e.g., see Goldman |2H]; the oldest 
reference appears to be from Hua |S5)- This allows us to compute in linear time a trapezoid r m 
that contains the vertical chord c m . 

Once T m has been identified, it is easy to compute x m (the x-coordinate of c x ). We desire the 
solution to the equation fiyy + q{x m ) = fJ*E + (A*( T m) — l{ x m)), where nw (resp., ^_e) is the area of 
all components of P \ r m that are adjacent to the left (resp., right) wall of r m , /i(T m ) is the area of 
r m , and q(x m ) is the area of the portion of T m to the left of coordinate x m . It is easy to see that 
q(-) is a quadratic function; thus, x m is readily computed as a root of a quadratic equation. Since 
fiw and fJ-E are readily computed in linear time once r m is identified, the computation of x m takes 
linear time in total. Similarly, we compute y m in linear time. □ 



6 Geodesic Li Distances in Polygons with Holes 

Now we discuss an even more complicated case, which arises when considering geodesic L\ distances 
in polygonal regions P that may have holes. Again, we analyze the set of locally optimal points: 
as long as a potential center can be moved in some axis-parallel fashion that lowers the average L\ 
geodesic distance to all the points, it cannot be optimal. 

The local optimality of a point Z is closely related to the subdivisions that it induces: for local 
optimality in the x-direction, the subdivision into W(Z) and E(Z) needs to be area-balanced; for 
local optimality in the y-direction, the subdivision into N{Z) and S(Z) needs to be area-balanced. 
(Refer to LemmaOl) The boundary between W(Z) and E(Z) is formed by bisectors in the shortest 
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Figure 14: A median trapezoid corresponds to a median in a node-weighted tree, with node weight 
representing trapezoid area. 



path map, SPM(Z) with respect to Z. It follows from basic properties of shortest path maps that 
the total complexity of this boundary is 0(n). (See, e.g., 07].) 

As we showed in LemmaHJ there is a neighborhood for each point Z € P in which the objective 
function / is cubic, provided that no bisector for Z meets a boundary vertex. This motivates the 
following lemma: 

Lemma 9 There is a subdivision of P of worst-case complexity I = 0(n 4 ), such that f is a cubic 
function within each face of the subdivision. 

Proof. Lemma0]implies that we are done if we can compute a subdivision of the claimed complexity 
such that we can move continuously between any two points in the interior of a connected cell of 
the subdivision, without any bisector encountering a vertex of the polygon during this motion. 
Provided that there is a position Z for which a bisector encounters a vertex v of the polygon, 
this vertex v is contained in W(Z) as well as in E{Z). Thus, there are two topologically different 
shortest paths from Z to v, one fully contained in W(Z), the other contained in E{Z). This 
implies that there are two topologically different shortest paths from v to Z, i.e., Z must lie on a 
watershed bisector in SPM(u). Therefore, the required subdivision is obtained by considering the 
0(n) watershed bisectors in each of the 0(n) shortest path maps with respect to polygon vertices. 
Each shortest path map has a complexity of 0(n), so the subdivision is defined by the overlay of 
0(n 2 ) line segments, yielding an arrangement of worst-case complexity / = 0(n 4 ). The example in 
Figure Hoi shows that even in the case of simple polygons, this bound on / is tight in the worst case. 
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(Chiang and Mitchell have studied similar arrangements that arise in overlaying shortest-path 
maps in the Euclidean shortest-path metric.) 



Q.(n) vertices 
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Figure 15: An example in which the overlay subdivision has complexity I = Q(n 4 



□ 

Considering the local optima for each cell of the arrangement allows us to obtain the following: 

Theorem 10 For geodesic L\ distances, a feasible point Z* = (x*,y*) in a polygonal region P 
with holes that minimizes the average (geodesic L\) distance f to all points in P can be found in 
worst-case time 0(1 + nlogn). 

Proof. The search for optimal solutions proceeds on a cell by cell basis, for each of the 0(1) cells 
in the overlay arrangement. The overlay arrangement can be computed in time 0(1 + nlogn), 
using known algorithms (|IJ E])- (We utilize the perturbation argument given in Section |2] in 
order to be able to assume, without loss of generality, that the bisectors are all polygonal curves, 
not regions of nonzero area.) For each cell, we spend constant time conputing the 0(1) candidate 
(local) optima inside and on the boundary of the cell. The function parameters for / within each 
cell can be determined in total time 0(1 + nlogn) by traversing the arrangement (e.g., by depth- 
first search in the planar dual graph of the arrangement) and doing simple 0(l)-time updates 
when changing from one cell to a neighboring cell. After determining the 0(1) candidate locations, 
we can determine a best among them by computing their objective values, again in total time 
0(I+n log n), by performing incremental updates to the objective function values during a traversal 
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of the arrangement. For any given cell of the arrangement, if there is a local minimum interior to 
the cell, the gradient V/ must vanish. Because / is the sum of two cubic functions, fi(x) and fe(y), 
within the cell, this means that we get a system of two quadratic equations (both components of 
the gradient must be zero) with two variables (x and y). Such a system can be solved in constant 
time using radicals. 

Similarly, we can determine the local optima with respect to variation along a boundary segment 
of a cell. For each segment, the gradient needs to be orthogonal to the segment. As in the straight- 
line case, this yields a quadratic equation that can be solved in constant time. 

Finally, there are 0(1) vertices in the arrangement, each of which we consider to be candidates. 

In total, then, we have examined 0(1) candidate local minima, in time 0(1 + nlogra). □ 

7 Multiple Centers 

We now discuss the fc-median problem of placing k centers into a polygonal region P, such that 
the overall average distance of all points p G P to their respective closest centers is minimized. We 
consider k to be part of the input and potentially large. 

Theorem 11 For polygons P with holes, it is NP-hard to determine a set of N centers that min- 
imizes the average geodesic L\ distance from the points in P to the nearest center. 

Proof. Our construction uses a reduction from Planar 3Sat, which was shown to be NP-complete 
by Lichtenstein [10]. We recall that a 3SAT instance I is said to be an instance of Planar 3SAT, 
if the following bipartite graph Gj is planar: each variable Xi and each clause Cj in / is represented 
by a vertex in Gj; two vertices are connected if and only if one of them represents a variable that 
appears in the clause that is represented by the other vertex. See Figure ITol The variable-clause 
incidence graph can be embedded in the plane without any crossing edges. 

First, the planar graph Gj corresponding to an instance I of Planar 3Sat with n variables 
and m = 0(n) clauses is represented in the plane as a planar rectilinear layout, with each vertex 
corresponding to a horizontal line segment, and each edge corresponding to a vertical line segment 
that intersects precisely the line segments corresponding to the two incident vertices. There are 
well-known algorithms (e.g., |54| ) that can achieve such a layout in 0(n) time and 0(n) space. See 
Figure We assume that the eventual layout is scaled appropriately by a factor of size G(n), 
such that the overall size is G(n 2 ). 

Next, the layout is modified such that the line segments corresponding to a vertex and all edges 
incident to it are replaced by a loop - see Figure El (top). At each vertex corresponding to a clause, 
three of these loops (corresponding to the respective literals) meet. Finally, the edges of any loop 
% are replaced by a sequence of 3q small squares (say, of size e = 0(l/n)) that are spaced apart at 
a constant distance (say, d = 0(1)) and are interconnected by narrow corridors (say, of width e 11 ) 
that have small enough total area that they do not greatly influence the overall average distance: 
The total area of all corridors is 0(n 3 *e n ) = 0(e 8 ), and the maximum distance between two points 
in corridors is 0(n 3 ), so the integral of pairwise distances over all corridor points is 0(e 5 ). Along 
each variable loop, the sequence of 3q squares is labeled "false" (index mod 3 in the sequence), 
"true" (index 1 mod 3 in the sequence), and "nil" (index 2 mod 3 in the sequence), in succession. 

Similarly, each vertex for a clause is replaced by a single small square and linked to the adjacent 
variable loops by three narrow corridors of length d, with adjacency encoding to the corresponding 
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Figure 16: The graph Gi for the Planar 3SAT instance / = {x\\lx2 VX3) N{x\ VX3VX4) A(x2 VX3 VX4), 
and its geometric representation. 



literal, i.e., connecting the clause square to a "true" square for an unnegated and to a "false" 
square for a negated literal. Note that no clause square is adjacent to a "nil" square. See Figure ITT1 
(bottom) for the overall picture. 

Let N = 3c = 0(n 2 ) be the total number of squares in all variable loops, and consider the 
placement of N centers. Because the total area of the corridors is small enough, neglecting them 
in the following discussion changes the resulting overall average distances by not more than 0(e 5 ). 
Furthermore, we assume without loss of generality that the placement of centers is locally optimal, 
so any center is placed as a median of the squares closest to it. It is readily checked (making use 
of the constant distance between adjacent squares and the negligible area of corridors) that this 
allows us to assume that all centers have been placed inside of squares. 

Now it is easy to estimate the overall average distance: We get an average distance of 



where n s d is the distance of the midpoint of square s to the midpoint of the closest square containing 
a center point; by construction, each n s is a nonnegative integer. As there can be at most c squares 
with n s = 0, we conclude that Yls=i n « — ^c, w ith equality if and only if each square is either 
occupied by or next to a square with a center. For this reason, we call a square s covered, if and 
only if n s < 1. It follows from the above description that there is a distribution of centers with 
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Figure 17: Replacing variables by loops (top); final polygon (bottom). Numbers 0, 1, 2 indicate 
the "true", "false" and "nil" squares. 



an average distance of ^ (@(e 4 ) + 2c) , if there is a covering positioning; on the other hand, we see 
that the average distance must be at least 2c ^ 1 if there is no covering positioning of centers. 

To establish the claim of NP-hardness, we show in the following that there is a satisfying truth 
assignment of the instance /, if and only if there is a distribution of centers to squares, such that 
all squares are covered. 

First assume that there is a covering. Consider the set of Cj "nil" squares in a variable loop. 
Clearly, no two of them can be covered by the same center; by construction, no "nil" square can be 
covered by a center in a clause square. This implies that each variable i requires precisely q centers 
to be placed in its squares, so all c centers must be placed on variable squares. As all 3cj squares 
in loop % must be covered by its Cj centers, and no center can cover more than three of the squares, 
we conclude that each center covers precisely three of the variable squares. Thus, all centers in 
a variable loop must be uniformally chosen to be all "true", all "false" or all "nil". Finally, each 
clause square must be covered from one of its adjacent variable squares, which means that setting 
variable i to the value indicated by the respective truth value satisfies the clause. Thus, all clause 
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squares can only be satisfied if there is an overall satisfying truth assignment /. 

Conversely, it is clear that for a satisfying truth assignment for instance /, placing Cj centers in 
the Ci "true" or "false" squares of variable i (corresponding to the truth setting of variable i) yields 
a covering distribution of centers. 

This concludes the proof. □ 

Note that the above proof can also be applied to the case of geodesic L2 distances, or when 
minimizing the maximum distance instead of the average distance; furthermore, the underlying 
proof technique can also be applied to other types of location problems. For example, see [21] f° r 
a game-theoretic scenario in which two players try to claim as much area as possible by placing 
centers, and the second player must place all of his points after the first player has played all of her 
points. 



8 Conclusion 

In this paper, we have given the first exact algorithmic results for the Fermat- Weber problem 
for a continuous set of demand locations. We have shown that for L\ distances in the plane, we 
can determine an optimum center in polynomial time, with the complexity ranging from 0(n) for 
the case of geodesic L\ distances in simple polygons, to 0(n 2 ) for straight-line distances in general 
polygonal regions, and 0(n ) for geodesic L\ distances in polygons with holes. Our results rely on a 
careful understanding of the local optimality criteria, together with the structure and combinatorics 
of shortest path maps. 



Extensions and Open Problems. 

(1) Our results can be extended to "fixed orientation metrics" defined by any constant number of 
directions. (The L\ metric is the special case in which the two fixed orientations are horizontal and 
vertical.) The local optimality conditions become more complex; however, the inherent algebraic 
complexity remains the same, for any metric whose disks are convex polygons. This extension 
allows one to approximate the Euclidean (L2) case to any desired degree of precision. 

(2) Our local optimality conditions generalize to the case of more general (non-uniform) nonneg- 
ative demand densities 5(p) by using the following observation. Regardless of the demand density 
function, any center location Z E P induces a subdivision of P into E{Z) and W(Z), and into N(Z) 
and S{Z) by shortest-path bisectors. Then the local optimality condition on Z requires that E{Z) 
and W(Z), and S(Z) and N(Z) are balanced in the following sense: instead of requiring that E{Z) 
and W(Z), as well as N(Z) and S(Z), have the same area, the balance condition is that locally op- 
timal points Z must have the integrals f pE yyr Z ) $(p)dp and f peE r z -\ 5(p)dp, and f pEN r z \ ${p)dp and 
IpeS(z) ${p)dp the same. Points with these properties are called 5-medians. Similar ideas can be 
used for describing boundary points. If, for a particular 5, there is a limited number of ^-medians, 
they can be computed in polynomial time, and it is possible to compare objective values in polyno- 
mial time, then we can determine a <5-center for the given region. This includes the case in which 
the demand function is given by point weights in combination with a uniform demand distribution 
over P, which is a problem formulated by Wesolowsky and Love |62j . It is also easy to see that 
the above methods can be applied for the case in which F ^ D and distances are straight-line 
L\ distances. Note that geodesic distances are not well-defined in this case; however, if we use a 
combination of straight-line distances outside of F and geodesic distances inside of F, our methods 
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still apply. 

(3) Our methods can also be applied in higher dimensions, by generalizing the local optimality 
conditions and carrying through the analysis in a very similar manner to the two-dimensional case. 
Figure ED shows that a generalization of Theorem however, does not hold in three-dimensional 
space (since any axis-parallel plane cuts the region into not more than two pieces), so we cannot 
use the same idea that allowed us in the two-dimensional case to exploit simplicity in achieving a 
better complexity than in the case of a polygon with holes. However, we can apply the technique 
of decomposing space into cells and studying the objective function within each cell. As in the 
two-dimensional case, the objective function is cubic for each coordinate, if P is a polyhedral 
region. 




Figure 18: In three-dimensional space, there may not be a feasible point that is a median of P in 
all coordinates. In this example, the intersection of the three median planes (orthogonal to the x-, 
y-, and z-ax.es) is a point approximately at the center of the bounding box, not lying within the 
solid. 

(4) Our methods for searching for local optima should be extendable to the case in which we 
have a constant (A; = O(l)) number of centers, e.g., k = 2. The centers induce a subdivision of 
P into several "Voronoi regions", corresponding to the set of points closest to each center. Each 
center must be placed optimally with respect to its region, which can be done by our methods. 
Thus, we are done if we have a suitable way to characterize the boundaries of Voronoi regions, 
which consist of a number of bisectors. While there are considerable technical details to establish, 
we believe that this approach will allow our results to generalize to multiple centers, and will lead 
to an algorithm of complexity 0(n ck ), for a small constant c. 

(5) It would be most interesting to discover an algorithm with worst-case complexity better 
than 0(n 4 ) that can compute an optimal center in the geodesic L\ distance for polygons with 
holes. Can we use geometric special structure to avoid examining all potential local minima in the 
cells of the overlay arrangement? 
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